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1. Introduction 

The bus transportation system in Cuernevaca, Mexico, has certain distinguished, innovative features and 
has been the subject of an intriguing, recent study by M. Krbalek and P. Seba in The situation is as 
follows. We quote from jT3j: In Cuernevaca there... 

is no covering company responsible for organizing the city transport. Consequently, con- 
straints such as a time table that represents external influence on the transport do not exist. 
Moreover, each bus is the property of the driver. The drivers try to maximize their income 
and hence the number of passengers they transport. This leads to competition among the 
drivers and to their mutual interaction. It is known that without additive interaction the 
probability distribution of the distances between subsequent buses is close to the Poisonian 
distribution and can be desribed by the standard bus route model .... A Poisson-like dis- 
tribution implies, however, that the probability of close encounters of two buses is high (bus 
clustering) which is in conflict with the effort of the drivers to maximize the number of 
transported passengers and accordingly to maximize the distance to the preceding bus. In 
order to avoid the unpleasant clustering effect the bus drivers in Cuernevaca engage people 
who record the arrival times of buses at signiflcant places. Arriving at a checkpoint, the 
driver receives the information of when the previous bus passed that place. Knowing the 
time interval to the preceding bus the driver tries to optimize the distance to it by either 
slowing down or speeding up. In such a way the obtained information leads to a direct 
interaction between buses... 
In ^31 J Krbalek and Seba describe their work in analyzing the statistics of bus arrivals on Line 4 close 
to the city center. They study, in particular, the bus spacing distribution and also the bus number variance 
measuring the fluctuations of the total number of buses arriving at a fixed location during a time interval T. 
Quite remarkably, Krbalek and Seba find that these two statistics are well modeled by the Gaussian Unitary 
Ensemble (QUE) of random matrix theory (RMT)(see figures 2 and 3 in ^Sj)- Our goal in this paper is to 
provide a plausible explanation of these observations, and to this end we introduce a microscopic model for 
the bus line that leads simply and directly to CUE. 

As noted in I13j the number variance for the buses is in good agreement with the CUE formula up to a 
time interval T=3 (see figure 3 in jl3|). As explained in |13j . this behavior is consonant with the fact that 
each bus driver, using the information given by the recorder, interacts with the bus immediately behind him 
and the bus immediately in front of him. In other words, the primary interaction is a three body interaction. 
As is well known, particle systems modeled by CUE involve interactions between all the particles. This 
means, in particular, that if we use CUE to model a system with nearest neighbor (or more generally, short 
range) interactions we should restrict our attention to statistics that involve only nearest neighbor (or short 
range) interactions. This is the case for the spacing distributions and hence one is able to account for the 
good agreement between CUE and the observed data for this statistic across the entire parameter range 
m PI : this is in contrast to the number variance where, as noted above, there is good agreement only up to 
T=3. 

The paper is organized as follows. In Section El we introduce our model. In Section El we describe the 
double scaling limits of interest and indicate how to analyze these limits for our model using standard 
asymptotic techniques from random matrix theory. Our main results on the spacing distributions, as well as 
on the number variance, are stated at the end of Section|31. In Section^lwe show how to modify the model 
to include alternative bus schedules. 
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Figure 1 . Conditioned Poisson walks arriving at location 3 



Our analysis is based on the fact that certain nonintersecting paths models lead to random matrix type 
ensembles, more exactly, orthogonal polynomial ensembles. Various manifestations of this phenomenon can 
be observed in [iiniiniEI- 
Finally, in view of the calculations that follow, we note that formula (4) in 13 for the density of the 
spacing distribution is the so-called "Wigner surmise" rather than the Gaudin distribution. The Gaudin 
distribution (see H28() below) gives the exact formula for the spacing distribution: However, the Wigner 
surmise is known to approximate this exact formula to high accuracy (see Mehta |15|^. 

Acknowledgments: The authors would like to thank Y. Avron for bringing the work of Krbalek and Seba to 
their attention. The authors would also like to thank P. Seba for providing useful information. The work of 
the authors was supported in part by NSF grants DMS0457335, DMS0402047, DMS0296084/0500923, and 
DMS0202530/05534 03, respectively. 



2. Basic Model 

The dynamics of our bus model takes place on the lattice Z as follows. Assume that there are n buses. 
Fix N > n and T > 0. At time t = the i*'' bus is at location 1 — i, 1 < i < n. The evolution of the bus 
process Lt — • • • , -^t"'), is given by n independent rate 1 Poisson processes conditioned not to intersect 

for t € [0, T], and subject to the terminal condition Lt = {N, N — 1, N — n + 1). Consider buses arriving 
at a fixed point x € {1, 2, TV — n + 1} with consecutive arrival times < ti < ■ ■ ■ < tk < T (see Figure 1). 

The probability density for the arrival times ti, t„, Ax{ti, tn)dti ■ ■ ■ dtn, is given by 

(1) 

P [{(0, 1 — i)} {(<i, X — 1)}; no intersection] P [{(ii, x)} {(T, N — {1 — i))}', no intersection] dti ■ ■ ■ dtn 

P [{(0, 1-i)} ^ {(T, TV - (1 - i))}; no intersection] ' 

The fact that the numerator in ^ splits into a product of two factors is a consequence of the strong Markov 
property of Lf. in particular, functionals of the parts of the paths {{ti,x) — > {T,N — (1 — j)), 1 < i < n} 
depend only on the increments of the process after some stopping times and hence, are independent of any 
functionals of the initial parts {(0, 1 — i) ^ {ti, a; — 1), 1 < z < n} of the paths. Note also that the endpoints 
{ti,x — 1), 1 < i < n, reflect the fact that the i*'' path jumps at ti, and not before. Similar dynamical models 
with T — oo were considered in |12j . 

As we now show, each term in the ratio has a determinantal form of Karlin-McGregor type ^Ol • 



Lemma 1. 



(2) P [{(0, 1 - i)} ^ {{ti,x - 1)}; no intersection] = det [ e *^ - — — 



{x + i~2)\ 

(3) P[{{ti,x)} {{T, N + (1 ~ i))}- no intersection] = det e^C^"*^'^- ^ rr 

\ {N + {I - i) — x)\ 



(4) 



' [{(0, 1-i)} ^ {{T, N+{1- i))}- no intersection] = det ( e"'^ 



{N + i-j)\ 



Proof. Since the proofs oflO and ^ equalities are similar, we only prove The proof is a variation on 
the standard Karlin-McGregor argument HP . Fix < ti < • • • < t„ < T and define the following stopping 
times: 



T2 



[m/{5e[0,ii]:LW=Lf }]Ar 
[m/{5G[0,i2]:4"=4''}]AT 



Tn-i = [zn/{,5G [0,i„-i] =Li")}] AT, 



Observe that 



det e"*^ 



(a; + i-2)! 



sg?i((7)]E 
E sgn{a)¥. 



llriA---Ar„_i=T + f-ri A---AT„_i=Tj<r 

J = l 



n-l 



(11, 



Ti A---Ar„_i = 



i=-.<r)n^{Lr.«=.-i} 



' [{(0, 1 — i)} — > {(ii, a; — 1)}; no intersection] 



(5) 



5] 5] sgn(a)E 



(llriA-Ar„_i=r,<T) ^{Lr.(''=x-1 



1} 



since the only nonzero term in the first sum is the term corresponding to the identity permutation. We need 
only show that the second summation vanishes. Let pj G S'„ be the transposition (j, j + 1). By the strong 
Markov property for Lt , 



sgn{a)^ 



(11, 



Ti A ■ ■ ■ A r„ 



.i=r,<T)n^{L-(')=x-l} 



+ sgn{pja)E, 



(11, 



Ti A---At„ 



2 — 1 



-1} 



0. 



Since for each j the action of the given transposition pj is an involution on S„ 



(llriA-..Ar„_i=Tj<T) Y\_ ^ 



{L,7'=x-1} 



This immediately implies that the second sum in equation vanishes leaving only the desired probability. 
The proof of the third equality (0)) is just the standard Karlin-McGregor argument JOI' In the proof of 
simply replace the r^'s by 



(6) n=TA inf{s e [0, T] : LI = 1'+^}, for i = 1, n - 1, 

a; — 1 by + (1 — i), and set tj = T for 1 < z < n. This completes the proof of the Lemma ^ 



□ 



The determinants in ^ and ||2J) are multiples of Vandcrmondc determinants. Indeed, simple computations 
show that 

Inserting ijTJl, and JHJl into J^l and setting yj = one obtains 
A{ti,...,tn)dti---dtn = A{yi,...,yn)dyi - ■ - dyn 

n 

where CN,n,x = (llILi (^+»-2)!(Af^-a:+(i-»))! ) ('^*'K(5v+bi)T)^^._^^ ■ The simple variable change yj = 
2yj — 1 induces the standard Jacobi weights on [—1, 1]: 

TL 

(10) CM,n,. Uiy, ~ y^f + - y,f'^-^"-'''dyi ■ ■ ■ dyn, 

where CN,n,x = CN,n,x- It is at this point that contact is made with RMT: If we view the y/s 

as eigenvalues of an n x n Hermitian matrix, then H10(l is precisely the joint probability density for the 
eigenvalues of matrices in the so-called Jacobi Ensemble . More precisely, for n x n Hermitian matrices 
{M = M*}, set 

(11) u;(Af) =det(l + Af)^-Met(l-M)^-^-"+Met(x(-.i,i)(M)). 

Note that w{M) is invariant under unitary conjugation of M, M UMU*, for all unitary U. Then 
Pn{M)dM = -^w{M)dM defines the Jacobi Unitary Ensemble , where dM = nr=i dM,i H^^. dM^j H-^j dM^^ 
and Zn is a normalization constant. A standard calculation (e.g. Mehta ^H]) then shows that the distribu- 
tion function for the eigenvalues, {A^}, of matrices in the ensemble is given by 1)101) with the identification 
yi — Xi- In summary, we see from pU|l that for finite N,n the arrival times for the buses at a location x 
are distributed like the eigenvalues of a random matrix from the Jacobi Ensemble. By techniques which are 
now routine in random matrix theory (see |15)). (|10|) is amenable to asymptotic analysis in certain double 
scaling limits. These limits (see Section O correspond to GUE: The fact that Jacobi Unitary Ensemble — > 
GUE in the double scaling limit is a particular example of the well known phenomenon of universality in 
RMT (see 1301111). 

In the previous calculations we considered the arrival times ti, ...,t„ for the buses at a fixed location, x. 
It is also of interest to consider the locations xi, at a fixed recording time t. It turns out that the 

distribution of the Xi^s is again given by a random matrix type ensemble, the Krawtchouk Unitary Ensemble. 
In Konig, O'Connell, and Roch discovered a similar relation between the Krawtchouk Unitary Ensemble 
and Poisson random walks conditioned not to intersect for all positive times. In analogy to Q), we have that 
the probability Ct{xi, ...,x„) for the positions xi > • • • > a:;„ of the buses at time t G (0,T) is given by 

P ({(0) ^ ~ i) ^ {t,Xi)}; no intersection) P {{{t, Xi) {T, N + 1 - i)}; no intersection) 
^ ' ¥ {{{0,1 - i) (T,N + no intersection) ' 

and by a Karlin-McGregor argument, similar to and even simpler than that in Lemma ^ we obtain the 
following: 

Lemma 2. 



/ ^Xj+i-l \ 

(13) P ({(0, l-i) ^ {t, X,)}; no intersection) = det e"*- -- 

V {Xj +1-1)1 J 



(14) 



' {{{t, Xi) ^ {T,N + 1 - i)}- no intersection) = det ( e^^"*^ 



(T-t) 



N+l-i-Xj 



(iV + l-z-x,)!/,.,-^i 



(15) 



P ({(0, 1 - i) -> (T, TV + 1 - i)}; no intersection) = det e 



Inserting ifT^ . and ((T5|l into ifT^ we are led to the following formula for £t(xi, x„) [5]: 

{N + i-iy.' 



Ct{xi, ...,x„) 



(16) 



where 



(17) 



[{N + n- 1)1] 

.1=1 

{N + n-l)\ 



n 



X 



n y^.!(Ar + „_i_y^.) 



(.-1). 

pW(l -p)^+"-l-!'j-, 



?i — 1, p = 4, and + — l>yi>---j/„>0. We have used the elementary identity 



det 



1 



n 



{N + i- 1)! 



to derive As indicated above, (|16(l corresponds to the Krawtchouk Unitary Ensemble. 

An essential aspect in the analysis of the random particle ensembles (I1U|I and (|16|) is that they are 
determinantah i.e. for each r = 1,2, the r*^ correlation function for the ensemble can be expressed in 
determinantal form det{K{xi, Xj))l for some appropriate correlation kernel K{x,y) (see for example, 
[151 1181 S El)- It turns out that the fixed time distribution l|16|l . in particular, can be extended to a 
dynamical random particle system which is also determinantal. Indeed, fix k times < ii < • • • < < T 
and let x'^) = {x[^\ ..^xi^^ }, I < j < k, denote the positions of the buses at times tj, 1 < j < k. Let 
Pj — Y 1 ^ 1^ j ^ ^- Then, using the Markov property of the bus system as above, we arrive at the following 
distribution for x^^\ ...,x'^^^: 



Prob[x 



(1) 



(fe)l 



const ■ Wi- 



,(1) 



i<3 



X det(F^,_p,(xf \xf)))det(Fp3_p,(a;f\xf )) • • •det(P^,,_,,_,(xf "^\a:f))) 



(18) 



{i-Pk) 



=-^1 (iV + n- 1 - xf )! 



where V^(x, y) = ^y^^y if 2/ > and Vp{x, y) ~ otherwise. Using a variant of the Eynard-Mehta theorem |S] 
as described in 3 , it follows directly that (|18|) is determinantal with an appropriate correlation kernel 
Kt-^tj{x,y). Further analysis yields the following block integral representation for the kernel: for x £ 



(19) 



Kti,tAx,y) 



(27ri)2 



{pj - (1 - Pj)t)y[l + t)^+n-i-yt-n ^s^t 
(k-(1-P,)s)^+K1 + s)^+' 



n—x o— n 



t 



where the integration contours are chosen as follows: s runs along a simple positively oriented contour which 
goes around and does not contain —1; t runs along a simple positively oriented contour which goes 
around 0; for pi > pj the s-contour contains the t-contour and for pi < pj the s-contour lies inside the 
f-contour. 

In the following section we will analyze (|10() and (|16() asymptotically as N,n, x,t,T oo. We plan to 
present the derivation and asymptotic analysis of (|18|l in a later publication. 
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3. The Double Scaling Limits 



We will first describe the arrival statistics of the buses (|lUf) in the limit N,n,x — ^ oo where ^ i/, 
— > 77, and < 1^,77 < 1. This scaling limit is the natural one from both the physical and mathematical 
point of view: the number of buses, n, should be proportional to the length of the bus route, N -\-n, and the 
arrival point, x, should not be too close to the beginning or the end of the route. The main result is that 
in the "unfolded" scale (see below), the Christoffel-Darboux kernel (equivalently, the two point correlation 
function) KN,n,x below, converges as — > 00, — > z/, — > 77 to the so-called sine kernel (see below), which 
is universal in random matrix theory. The proof of the convergence KN,n,x — > ICoo for the Jacobi ensemble 
was given first by Nagao and Wadati 'IB'. 

In RMT the analysis of the spacing distribution proceeds via the analysis of the gap distribution (see, for 
example, [T^f7]l: 

(20) Fn.u ((c, d)) = Prob(yj i (c, d), i=l,...,n). 

There are three steps involved in the analysis of Pn.u ((c, 0?)): (a) The first step is to express VN,n{{c, d)) in 
Fredholm determinantal form, 

(21) P^,„((c,d)) =det(/-if)i2[,,<„, 

for some kernel operator K expressed in terms of certain appropriate orthogonal polynomials (see, in par- 
ticular, Tracy and Widom ^H]); (b) the second step of the analysis is to determine the appropriate scalings 
for c and d; this is accomplished by evaluating the so-called equilibrium measure for the associated loga- 
rithmic potential theoretic problem ^| Ej ; (c) the third step involves a detailed asymptotic analysis of the 
appropriate associated orthogonal polynomials in the double scaling limit. 

(a) For j = 1, 2, ... , let pf'"'^(y) — 'jf'^'^y^ + 7^^'"'"^ > denote the Jacobi polynomials obtained by 
orthonormalizing {1, y, y^, ...} with respect to the weight WAr,„,x(y) = (1 — y)^~("+^~i) (1 -|- ?/)^~^X[-i.i] (?/)• 

Set (j)f'^'^{y) = P^'^"^{y)wjs! n xiy)^ ^'^^ j ~ 2: Then, the probability that there is no bus arrival in the 
interval of time [aN,n,x,bN,n,x], ^N,n{{0'N is given by 

(22) det{I - ifjV,n,a;)L2([ajv,„,x,&N,„,x])' 

where 

(23) Ajv,„,a;(z, Z ) = In ; , 



and l„ is given by the formula -j^^ y %^ where hn and fc„ are constants arising in the three term recurrence 
relation for general orthogonal polynomials; for the Jacobi polynomials, hn and fc„ are given, respectively, 
by formulas 22.2.1 and 22.3.1 of Abramowitz and Stegun in the case a = N — {n + x — 1) and f3 = x — 1. 

(b) By standard methods (see Saff and Totik |17| and Deift 0), the equilibrium measure can be determined 
explicitly for the problem at hand. It turns out that the support of the equilibrium measure is an interval 
[a, b] C [—1,1], and for fixed i/, 77 > satisfying v + rj < 1, the measure takes the form: 



, / X , V\/ (x — a)(6 — x) 

24 i:(x)dx = , ^ ' = dx, 

7rV(l + a)(l + 6)(l-a;2) 

where a and b satisfy the relations: 

V _ (7/ + 7^-1) , _ V (^ + ^^-1) 



V(l + a)(l + fo) ^{l-a){l-by ^{l + a){l + b) ^{l-a){l-by 

and the radicals represent the positive square root. In the symmetric case where 77 — ^^j^. 




b=-a= 

(c) Let To G (a,fe). A standard calculation in RMT show that the expected number of particles per unit 
interval in a neighborhood of tq « mpiro). Thus, changing scales f 1— > s := nip{To)t, we see that the expected 
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number of particles per unit s-interval is 1: this process of rescaling is known as "unfolding" the data. The 
probability that there are no bus arrivals in [tq — ti,{To)n ' 4>{To)n '\ becomes 

(25) rfet(/-/C„)i2(_,^,), 
where 

,,,, r r . ^ / '^^ ^""("° + Jik^^^'^n-riro + - 0,T-"'-(ro + ^)e""(ro + ^) 

(2b) AL„(To,4,pj = /„ . 

Analyzing the asymptotics of the Jacobi polynomial as iV, n, x —> oo as above, we finally see that 
(27) det (/ - /C„)^2(_,^^) det {I ~ /Coo)i2(_,^^) , 

where /Coo(C,p) = ^'^^ so-called sine kernel. To compute the limiting conditional probability that 

given a bus arrival at time s, the next bus arrives at time t we simply compute (see, e.g. |19p 

(2^) i^^oo W)^t'^"' ^^)lHs,) - ^^^^^ - ^-)^^(.*) • 

Formula (|28|l . the Gaudin distribution, gives the exact spacing distribution of RMT; as noted in the Intro- 
duction, this distribution in known to be well approximated by the Wigner surmise used by Krbalek and 
Seba |T5j to analyze their observations of the bus arrival times. 

In addition, using the fact that K]\[^„,x — > in the unfolded scale, one easily sees that the number 
variance, HN^n,xi'), say, converges i.e., 

(29) H^{s) := lim Hn.u.x ( — 

exists. Moveover, Hao{s) is precisely the number variance for GUE, as given, for example, in formula (16.1.3) 
in iiSj. As s becomes large, Hr^{s) ~ ■^^(log(27rs) -t- 7 -I- 1) where 7 is Euler's constant, as noted in formula 
(5)of[TS|. 

The analysis of the fixed time bus locations (|16|) as TV, T, i ^ 00, y ^ Ci, ^ C2, and ^ C3 is very 
similar and follows the same general procedure (a), (b), (c) above. We note that the Krawtchouk ensemble 
was first analyzed in the same scaling limit by Johansson in his analysis of the Aztec diamond |^. Again, 
one finds that in the limit the statistics of the bus locations is governed by the sine kernel. 

4. The Circular Bus Route 

In the previous bus model, we considered the buses moving from one terminus at i = to a second 
terminus at time t = T. In this final section, we consider a bus model on a circular route. Konig, O'Connell, 
and Roch(see again |^) investigated a related model on the discrete circle. Let Zm be the discrete circle 
with nodes labeled by {0, 1,2, ...,M — 1}. We will consider the case of fc < Af buses traveling along Zm- 
The buses start at positions < iVi(O) < • ■ • < Nk{0) < M — 1 and evolve as independent Poisson processes 
conditioned not to intersect. Let 9,9 £ Zj\,f and note that the transition probability on Zjv/ for a single rate 
1 Poisson process to travel from state 6* to is 



(30) Pt{9,9)^e 

I— — CO 



y. ■ 

^ {9^9 + lM)\ 



Let {9i, 9k) be a fc-tuple of distinct elements of Zm such that < 6*1 < ... < ^fc < M — 1. Let (6*1, 9k) 
be another fc-tuple of distinct elements of Zm such that there exists a cyclic permutation a £ Sk for which 
(0^(1), 6'(T(fc)) satisfies < 9^(i) < ■■■ < 9^(k) < M—1. Given that the buses begin at positions 0i < ... < 9k 
at time 0, the probability that they are at locations 9i, ...,9k at time t, and have not intersected in the mean 
time, is given by the determinantal expression 

(31) det (pt{9,,9^)Y . 

We prove this expression by adapting an alternate proof of the Karlin-McGregor formula from the line to 
the circle as follows. 
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Let be defined as the subset of Jj^ /MlJ^ for which there exists a cychc permutation a £ Sd such 

that < < •Z(t(2) <i • • • < < A/ — 1. For / : ifl.jy M, define the operator 

k 

(32) i/(z) = ^[/(zi, + 1, zfe) - /(z)] , for z e Z^,,> - aZ^,,>. 

As we will see, in order to prove the circular Karlin-McGregor formula H31|) , it is sufficient to show that the 
function gg : [0,t] x Z^j> R defined by gg{t,6) = det{pt-s{di,9j)) is the unique solution to the equation 

(33) (ds + L)gg = 0, 

(34) 5el[o,t]x9z^j> = 0, 

(35) gi^{t,e) = Sg{e). 

One checks easily that gg is a solution of H33|) . In order to check that the solution is unique, one needs to 
prove a maximum principle for this equation, but this is easily done by mimicking the standard proof for 
parabolic equations. 

We now show that gg{0,6) gives the desired probability. Let N{s) — {Ni{s), Nk{s)) be the Poisson 
process on Z|^> with initial condition A^(0) = (6*1, ...,6lfc)- Let r = inf{.s e [0,t] : N{s) e 9Z^^>}. By Ito's 
formula, 

(36) g^(^s,N{s/\T)-gg{O,0)- r\dr + L)gg{r, N{r))dr 

Jo 

is a martingale. Taking expectations, we obtain 

(37) gg{0, 6) = ¥Fg{s, N{s A r) for s G [0, t]. 

Letting s — t, 1)3 7|) becomes g§{0, 9) — K^g{t, N{t A r). Using the definition of r and the terminal condition 
3511. we see that 



(38) 5g(0, 61) = ¥'^(t < T and Nt = 0), 

which proves (|31|) . Thus, for buses starting at locations < < • ■ • < 6*^ < M — 1 at time t = 0, we have 
Prob(busses are at locationsO < 6*1 < • • • (?fe < Af — 1 at time t ; no intersection for all < s < t) 
= ¥{t<T and N{t) =9)= gg{0, 0) = det{pt{0^, 0,)%^^, 

as desired. 

This formula enable us to obtain, in particular, a Karlin-McGregor type formula for the solution of the 
following natural problem for the buses on a circular bus route. Imagine that the initial locations of the 
buses are &i 0i — i ^ \ for i = 1, fc. Suppose that the buses return to these locations at some fixed time T 
later. For any time < i < T, it immediately follows from the formula above and the Markov property that 
the distribution, Qt{0i, ■■■,0k), of the locations of the buses conditioned on arriving at 0i, ■■■■,0k at time T, 
is given by the formula: 



(39) Qt{0i,^^;0k) 



detipti0^,0J)det{pT-t{0^,0J)) 



det{pTi9^,9j) 

into which we may now substitute (|31|l . Observe that the number of rotations about the circle is not fixed. 
We plan to investigate the asymptotic behavior of this model in a future paper. 
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